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Abstract 



We calculate the perturbative value of the Free Energy in Lattice QCD in 
three dimensions, up to three loops. Our calculation is performed using the 
Wilson formulation for gluons in SU(N) gauge theories. 

The Free Energy is directly related to the average plaquette. To carry 
out the calculation, we compute the coefficients involved in the perturbative 
expansion of the Free Energy up to three loops, using an automated set of 
procedures developed by us in Mathematica. The dependence on N is shown 
explicitly in our results. 

For purposes of comparison, we also present the individual contributions 
from every diagram. These have been obtained by means of two independent 
calculations, in order to cross check our results. 

Keywords: Lattice QCD, Lattice perturbation theory, Free Energy. 
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I. INTRODUCTION 



In this paper we calculate the Free Energy of QCD on the lattice, up to three loops in 
perturbation theory. We are interested in the three dimensional case, involving only gluons. 

The choice of studying the Free Energy is motivated by the fact that, apart from be- 
ing a simple physical observable which can be calculated in perturbation theory, it is used 
to determine the static thermodynamic properties of a physical system. Furthermore, this 
quantity is ideal for studying various aspects of QCD. For example, the Free Energy is often 
employed in the context of thermal QCD, for the purpose of characterizing the deconfine- 
ment transition between the low temperature regime, ruled by confinement, and the high 
temperature regime ruled by asymptotic freedom. Asymptotic freedom guarantees a small 
coupling constant at large temperatures, which in turn allows for a perturbative evaluation 
of observables. In addition to that, the perturbation expansion of the Free Energy is involved 
in the determination of the "finite part" of the gluon condensate in lattice regularization. 

Finally, the plaquette expectation value, which is directly related to the Free Energy, is 
considered in the study of an effective 3-D pure gauge theory called Magnetostatic QCD 
(MQCD) [1], which is matched to QCD using perturbation theory. Due to the superrenor- 
malizable nature of the 3-D theory, only a finite number of divergences arise and these can be 
computed perturbatively, allowing for a complete match between the lattice regularization 
scheme and a continuum scheme such as MS. MQCD is an instance of "dimensional reduc- 
tion", whereby the static properties of a (3+l)-dimensional field theory at high temperature 
can be expressed in terms of an effective field theory in 3 space dimensions. Dimensional re- 
duction has been applied to QCD (see, e.g., [2,3]), in order to resolve a longstanding problem 
regarding the breakdown of the perturbation expansion for the Free Energy [4]. 

The evaluation of the Free Energy in four dimensions has already been carried out to 
2 loops [5] and 3 loops [6,7]. In three dimensions, the only results available so far have 
been obtained through the method of stochastic perturbation theory [8]. Our results can 
be used both for comparison with existing results extracted from different methods, or 
for calculations involving the perturbative expansion of the Free Energy. With respect to 
this, the coefficients of the perturbative expansion allow the subtraction of all the divergent 
contributions from the gluon condensate [9] and improve the previous estimates of Ref. [8]. 



II. CALCULATION 

In standard notation the Wilson action for gluons reads: 

Sl = X E Tr[l-U^(x)} (1) 

i/0 x, /u, v 

Here U^ v {x) is the usual product of link variables U^(x) along the perimeter of a pla- 
quette in the \x-v directions, originating at x; go denotes the bare coupling constant. Powers 
of the lattice spacing a have been omitted and may be directly reinserted by dimensional 
counting. 

We use the standard covariant gauge-fixing term [10]; in terms of the gauge field Q^{x) 
[U^(x) = exp^goQ^x))}, it reads: 
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Sgf = A E Tr {^QMAzQ v (x)}, A~Q v (x) = Q v (x - p) - Q„{x). 



(2) 



Having to compute a gauge invariant quantity, we chose to work in the Feynman gauge, 
Aq = 1. Covariant gauge fixing produces the following action for the ghost fields uj and uj 



S gh = 2^Tr{(A>(x)) t (A^(x)+ lff „ [Q,{x)^{x)\ + % -f [Q,(x) : A+u>(x) 

Q,{x)\Q^x),A+uo{x)]] 
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Qn(x), [Qn(x), [Qh(x), [Qh(x), A+u(x) 



+ 



A+u(x) = uj{x + p) — ui{x). 



(3) 



Finally the change of integration variables from links to vector fields yields a Jacobian 
that can be rewritten as the usual measure term S m in the action: 



S m = E{^| Tr {(g M (x)) 2 } + ^| Tr {(Q,(x)) 4 } + ( T r {(Q,(x)) 2 }) 



(4) 



In Sgh and S m we have written out only terms relevant to our computation. The full 
action is: 



S — SL + Sgf + Sgh + S m . 



(5) 



In D-dimensions, the average value of the action density, S/V, is directly related to the 
average lxl plaquette P: 



where: 



(s/v) = D{D n 1) He{p)). 



£(P) = l--Ue[Tr(P)], 0=™ 



We will calculate (E) in perturbation theory: 

(E) = ci gl + c 2 #o + c 3 gl + 



(6) 



(7) 



(8) 



The n-loop coefficients c n have been known for some time up to 3 loops in the 4- 
dimensional theory [6]. 

The calculation of c n proceeds most conveniently by computing first the Free Energy 
— (\nZ)/V, where Z is the full partition function: 



[VU]exp(-S). 



(9) 



The average of E is then extracted as follows: 
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(E) = ~ 



D(D-l)^- 1 d f\nZ" 



(10) 



In particular, the perturbative expansion of (lnZ)/V is: 

QnZ)/V = do - (iV 2 -l) ln/3 + | + | + 



(11) 



We can now write the perturbative expansion of (E), by combining Eqs. (10) and (11): 



(E) 



D(D-l) 



(D-l) (N 2 -l) ck 2d, 

2 (3 (3 2 (3 3 



The coefficients c n involved in Eq. (8) can be specified by the general expression 

_/ £>(£> (n-lR-x 
Cn "l 2 J (27V)n ' n ^ 

which leads immediately to the following relations for D = 3: 

c 2 = rf!/(12iV 2 ), 
c 3 = d 2 /(12N 3 ). 



(12) 



(13) 



(14) 
(15) 



III. RESULTS 

A Total of 36 Feynman diagrams contribute to the present calculation, up to three loops, 
as shown in Figure 1. The involved algebra of the lattice perturbation theory was carried 
out using our computer package in Mathematica. The value for each diagram is computed 
numerically for a sequence of finite lattice sizes L, followed by an extrapolation to L — > oo 
(see next Section). The finite- L results, on a per diagram basis, are available upon request 
from the authors. We performed two independent calculations for each diagram, in order to 
verify our results. 

Certain diagrams, considered individually, are infrared divergent; they contain subdia- 
grams which renormalize the gluon propagator to one loop, as shown in Figure 2. Taken 
as a group, these diagrams become infrared finite and their value can be extrapolated to 
infinite lattice size; extrapolation leads to a (small) systematic error, which is estimated 
quite accurately. 

The coefficients q involved in the perturbative expansion of (E) can be calculated by 
summing all contributing diagrams. They take the form: 



N 2 -l 






6N ' 








:) 




t 3iV 


1} 


, 7V2 + C 3,l + C 3,2 iV 



(16) 
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21 



22-36 



Fig. 1 Diagrams contributing to the Free Energy at one loop (1), two loops (2-6) and three loops (7-36). 
Solid (dashed) lines represent gluons (ghosts). The filled square is the contribution from the "measure" 
part of the action. Black bullets stand for an effective vertex, which is shown in Fig. 2. 
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Fig. 2 Effective vertex involved in diagrams 22-36. The filled square is the contribution from the 

"measure" part of the action. 



The index j runs over all contributing diagrams. The coefficients c^|, are pure 
numbers and we calculate them for each diagram separately. Their numerical values are 
listed in Tables I and II. For some diagrams, the coefficients are known analytically. The 
expressions are presented as a function of the 1-loop vacuum diagram Pi (see below): 

eg = (1 + 28Pi - 72P 2 )/48 

4 5 J = -A/12 

eg = -PJ8 

eg = (504Pi - 1296P 2 - 6)/5184 

cH = (3240Pf - 792P 2 - 63Pi - 1)/5184 

eg = -Pf/288 
c { S = -AV384 
4 1 ? = -P 2 /288 

(17) 

The value of Pi can be expressed in terms of the Bessel function of the first kind I : 



dt 



* dp -t4sin 2 (p/2) 
-* (2tt) 



o 



-2t 



J (2t) 



1 L (27T)3 3 J Q 

4]Tsin 2 (^/2) 

= 0.2527310098590 •• • (18) 



The ghost loop in diagram 21 leads to a vanishing contribution by color antisymmetry. 

Using our results we can express the average Free Energy as a function of go, for specific 
values of N: 

N = 2: (E) = (1/4) g 2 Q + 0.01453916(1) #£ + 0.0053459(2) #° + ••• (19) 
N = 3: (E) = (4/9) # 2 + 0.05420318(2) ^ + 0.0317648(7)^ + ••• (20) 



Presented in Fig. 3 is the behaviour of the average Free Energy as a function of 1/(3, for 
N = 2 and iV = 3. 



IV. DISCUSSION AND CONCLUSIONS 



Several remarks are in order regarding our computation. 

• As mentioned above, certain subsets of diagrams are infrared finite only when summed 
together, and they correspond to insertions of the one-loop renormalized gluon propagator. 
Similarly, diagrams 11, 12 and 13 renormalize the ghost propagator; however, they are 
separately finite. 

• The manipulation of the algebraic expressions in order to reach a code for numerical loop 
integration, is completely done in Mathematica. The production of vertices as well as their 
contraction is fully automated. A number of routines render automatic the simplification of 
the expressions (such as color contractions and Lorentz index manipulation) and exploitation 
of their symmetries. We have also developed an efficient three dimensional "integrator" to 
automatically convert the integrand to Fortran code. The size of a typical input integrand 
is up to < 1000 terms after simplification. Finally, the extrapolation of data to infinite L is 
also automated. The only part which is done by hand is the enumeration of diagrams, but 
this is simple enough to do, even up to 4 loops. Each n-vertex diagram is represented by 
an n x n "incidence" matrix M, where My is the number of gluons joining vertices i and 
j and by a similar matrix for ghosts. The procedure described above is not CPU intensive 
as regards symbolic manipulations in Mathematica. Of course, this is due to the usage of 
Wilson action for gluons (improved actions are more cumbersome in terms of CPU as well 
as RAM). 

• Numerical integrals are performed for finite lattice sizes 1 , up to L = 100 3 , and the results 
are subsequently extrapolated to L — > oo. Diagrams 18-20 (those having the form of the 
Mercedes emblem) are more CPU intensive and require approximately 40 hours each (for 
L = 40 3 ) on a recent PentiumlV computer (3.2 GHz); for these diagrams, integration over 
the three loop momenta can only be performed in a nested fashion, rather than in parallel. 
Fortunately, lattice sizes up to L = 40 3 already lead to a sufficiently small systematic 
error in these cases. For the rest of the diagrams, two out of the three loop momenta can 
be integrated in a parallel fashion, and the calculation requires approximately 1 hour per 
diagram for L = 100 3 . 

Extrapolation to infinite lattice size is of course a source of systematic error. To estimate 
this error, our procedure carries out the following steps: First, different extrapolations of our 
numerical data tl are performed using a broad spectrum of functional forms f k {L) (around 
50 of them) of the type: 

f k (L) = Y,et ) L-^\nLy (21) 

A total of N k coefficients appear in the k th functional form; these coefficients are deter- 
mined uniquely using the results on N k lattices of consecutive size L. 



X A detailed treatment of the zero modes on a finite lattice is not necessary for the purposes of 
extrapolation to L — > oo 
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For the k th such extrapolation, a deviation dk is calculated using alternative criteria for 
quality of fit. One possible criterion, as an example, is the difference: dk = \f k (L") — r L * |, 

* (k) 

where L is the largest lattice size which was not used in the determination of e\j . Finally, 
these deviations are used to assign weights d^ 2 / (Y,k d^ 2 ) to each extrapolation, producing 
a final value together with the error estimate. We can check the reliability of our error 
estimates in a number of ways: In those cases where the result for a diagram is also known 
analytically, this result coincides with the numerical one within the systematic error; also, 
extrapolations which incorporate new data from larger lattices are compatible with previous 
results, again within systematic error. 

The 3-dimensional case is slightly different and more delicate than the case of 4 dimen- 
sions; in the latter only even powers of 1/L must be used (i = even in Eq.(21)), along with 
possible logarithms. In 3 dimensions one must also use odd negative powers of L, since 
the main difference between a sum over discrete values of momenta and the corresponding 
integral comes from omitting the contribution of the zero mode. This quantity behaves like 
I = J v d 3 p/p 2 where V is a sphere of diameter 2n/L . We see that I = An 2 /L. 

Using larger lattices is of course computationally more feasible in 3 dimensions. At the 
same time, it is also more necessary than in 4-D, since the infrared behaviour of loop integrals 
worsens in 3-D. Suffice it to recall that 4-loop integrals in three dimensions will bring about 
a genuine divergence, requiring introduction of an IR regulator. 

• Referring to diagrams 22-36, integration in 3-D with use of an IR regulator (in our case, 
L) may give extra terms which do not exist in the 4-D case, leading to potentially wrong 
results, when the regulator is taken to its limit value (L — > oo). Due to the presence of a 
loop involving two propagators at the same momentum, each of these diagrams exhibits a 
behaviour of the type: 

(a+y)(cL + d) (22) 

After summing all diagrams and extrapolating to the limit L — > oo, the result we get 
has the form ad + be (terms like acL add up to zero). Of course, the terms be are spurious 
and must not be taken into account. A similar procedure in 4 dimensions leads instead to: 

(a + -^r)(c\nL + d) (23) 

so that spurious terms are absent: (b/L 2 )(c In L) — > 0. To solve this problem, we can 
write the mathematical expression for each vertex involved in diagrams 22-36 (Fig. 2) in a 
subtracted way: f(p) = /(0) + (f(p) — /(0)), where p is the momentum on external lines. 
Once all the subtracted vertices are summed, the terms /(0) cancel among themselves. In 
addition to that, the terms f(p) — /(0) are at least of first order in p and factors of the type 
cL involved in Eq. (22) cannot arise. The subtracted diagrams take the form (a + b/L)(d) 
and thus, extrapolation gives the desired results ad. 

• For purposes of comparison, we can express the perturbative expansion of the Free Energy 
shown in Eqs. (19) and (20), as a function of 1/(3: 

N = 2: (E) = (i) +0.2326265(2) ^ + 0.34214(1) (^) 3 + --- (24) 
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N = 3 : (E) = (8/3) (±) + 1.951315(1) (±) 2 + 6.8612(2) (I)' + • • • (25) 

For the case N = 3, the value (2iV) 3 c 3 = 6.8612(2) improves a previous estimate 
(2iV) 3 C3 = 6.90^g t j > 2 [8], which was obtained with the method of stochastic perturbation 
theory; it allows for a more stable extraction of the "finite" part of the gluon condensate, 
as expressed in Ref. [9] (using the notation of that reference for the symbols q): 

where the logarithmic coefficient c 4 is known from the continuum vacuum energy density in 
the MS scheme: 

* = gii 7#(§-^). ^=3) =2.92942132... (27, 
and C4 requires a 4-loop calculation in lattice perturbation theory. 

• A possible extension of the present calculation to 4 loops is feasible with the techniques 
used in this paper, but one must take into account some additional difficulties: First of all 
there is a total of 250 diagrams to be calculated (43 pure gluon diagrams, 115 diagrams 
with one ghost loop, 37 diagrams with two ghost loops, 4 diagrams with three ghost loops 
and 51 diagrams with measure terms) and computation time is proportional to L 12 (in most 
diagrams, loop integrals are factorizable, with the exception of a "nonplanar" diagram and its 
two ghost variants). Furthermore, due to the existence of a genuine logarithmic divergence, 
one must insert an IR cutoff m leading to results of the form cm(m) + d. Of course, the 
vertices and the contracted expressions, as well as a 4-loop version of our integrator, can 
still be produced automatically in this case. 

Acknowledgements: This work is supported in part by the Research Promotion Foun- 
dation of Cyprus (Proposal Nr: ENTAS/0504/11). 
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TABLES 



TABLE I. Pe 


r-diagram 


contributions to oi- 


Diagram 


Co n 


Co 1 


2 


-1/12 


0.0724503107... 


3 





0.03394382(1) 


4 





-0.00383019(1) 


5 





-0.0210609174... 


6 





-0.0315913762... 


Sum 


-1/12 


0.04991165(2) 





rp A T > T T- IT T ) 1* 

lAtJLlij 11. Fer-dia 


gram contributions to 


C3- 


)iagram 


C3,0 


c 3,l 


C3,2 


7 


-1/216 


0.00744542215... 


-0.0029334803.. 


8 








-0.0002217811.. 


9 





0.01131461(2) 


-0.01213765(2) 


10 








-0.0001663358.. 


11 








-0.0002217811.. 


12 








-0.00016133(1) 


13 








-0.00011990(1) 


14 








-0.00164058(1) 


15 


0.0094117909(2) 


-0.00313726(1) 


0.00306459(1) 


16 








-0.00011644(1) 


17 








-0.00008704(1) 


18 








0.00128600(6) 


19 








-0.00008896(2) 


20 








-0.00004025(1) 


21 











22-36 


1/72 


-0.03170827(1) 


0.01911231(5) 


Sum 


0.0186710502(2) 


-0.01608550(3) 


0.00552737(9) 
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